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ABSTRACT 


This  report  develops  by  matrix  methods  the  integral  solution 
of  the  wave  equation  for  point  sources  of  harmonic  waves  in  a 
liquid  layer  of  a  multilayered  liquid-solid  half  space  in  a  form 
convenient  for  numerical  computation  on  a  high-speed  digital 
computer.  Only  the  case  is  considered  of  a  high-speed  liquid 
bottom  underlying  the  stack  of  layers.  The  integral  over  wave 
number  has  singularities  in  the  integrand  and  is  conveniently 
transformed  into  the  complex  plane.  By  a  proper  choice  of 
contours,  complex  poles  are  displaced  to  an  unused  sheet  of  the 
two-leaved  Riemann  surface,  and  the  integral  solution  for  the 
multilayered  system  reduces  to  a  sum  of  normal  modes  plus  the 
sum  of  two  integrals,  one  along  the  real  axis  and  the  other 
along  the  imaginary  axis.  Both  integrals  are  evaluated  by  a 
Gaussian  quadrature  formula.  Sample  computations  are  presented 
for  low-frequency  propagation  in  the  Arctic  Ocean  sound  channel. 
These  are  preliminary  computations  and  the  ice  layer,  which 
averages  three  meters  in  thickness,  is  not  included  in  the 
layered  system.  The  effects  of  the  ice  layer  on  propagation 
are  currently  under  investigation. 


2. 


FIGURES 


Figure  1. 
Figure  2. 
Figure  3 . 

Figure  4. 

Figure  5. 


Figure  6. 


Figure  7. 

Figure  8. 
Figure  9. 
Figure  10. 


Multilayered  half-space. 

Contours  for  integration. 

Variation  of  sound  velocity  with  depth  for  Model  A. 

Table  2  gives  additional  parameters  for  this  model. 
Computations  for  Model  A.  Phase  -  and  group-velocity 
dispersion  and  excitation  function  of  pressure 
dependent  only  on  layering. 

Computations  for  Model  A.  Variation  with  range  of  the 
absolute  value  of  pressure  for  the  normal-mode  contri¬ 
bution  of  the  sound  field.  Source  depth  150  m.  Hydro¬ 
phone  depth  50  m.  Source  frequency  10  Hz.  Source 
pressure  amplitude  1  dyne/cm  re  lm. 

Computations  for  Model  A.  Variation  with  range  of 
absolute  value  of  pressure  for  the  integral  contribution 
of  the  sound  field.  Source  and  detector  same  as  for 
Figure  5. 

Computations  for  Model  A.  Variation  with  range  of  absolute 
value  of  pressure  of  the  total  sound  field.  Source  and 
detector  same  as  for  Figure  5. 

Same  as  Figure  5  but  computations  carried  to  longer  range. 

Same  as  Figure  6  but  computations  carried  to  longer  range. 

Same  as  Figure  7  but  computations  carried  to  longer  range. 
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DEFINITION  OF  SYMBOLS 

compressional-wave  velocity  in  the  m-th  layer 

shear-wave  velocity  in  the  m-th  layer 

thickness  of  the  m-th  layer 

vertical  coordinate 

range  between  source  and  detector 

time 

angular  frequency 
phase  velocity 
wave  number 

imaginary  part  of  a  complex  number  x  +  iy 
velocity  potential  in  the  m-th  layer 
normal  stress  in  the  m-th  layer  parallel  to  z  axis 
pressure  in  the  m-th  liquid  layer 

horizontal  particle  displacement  in  the  m-th  layer 
vertical  particle  displacement  in  the  m-th  layer 
horizontal  particle  velocity  in  the  m-th  layer 
vertical  particle  velocity  in  the  m-th  layer 
density  in  the  m-th  layer 
density  at  the  source. 


;  Bessel  function  of  order  0 

;  Y  bessel  function  of  order  0 

»  O 

;  K  Bessel  function  of  order  0 
;  Hankel  function  of  the  first  kind  of  order  0 


(A');  Hankel  function  of  the  second  kind  of  order  0 
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INTRODUCTION 


This  report  develops  the  integral  solution  of  the  wave 
equation  for  point  sources  of  harmonic  waves  in  a  liquid 
layer  of  a  multilayered  liquid-solid  half  space  in  a  form 
convenient  for  numerical  computation  on  a  high-speed  digital 
computer.  Only  the  case  of  a  liquid  half-space  underlying 
the  stack  of  layers  is  considered.  Furthermore,  it  is  assumed 
that  the  speed  of  sound  in  this  last  layer  of  infinite  thick¬ 
ness  is  greater  than  the  speed  of  sound  in  liquid  layers 
immediately  above  it.  The  problem  of  immediate  concern  is 
low-frequency  propagation  from  harmonic  sources  in  the  central 
Arctic  Ocean  at  ranges  up  to  ten  times  the  water  depth  in 
shallow  and  intermediate  water  depths.  Solutions  are  develop¬ 
ed  for  the  pressure  perturbations  detected  by  a  hydrophone  at 
depth  or  the  ice  vibrations  detected  by  a  seismometer  on  the 
ice  surface.  A  modification  of  the  formulation  is  possible  to 
include  harmonic  vibration  sources  on  or  in  the  ice,  although 
it  appears  at  the  present  time  that  surface  sources  are  of 
limited  application  (Greene,  1968). 

The  solution  of  the  wave  equation  presented  here,  based 
on  the  Thomson-Haskell  matrix  method  (Thomson,  1950;  Haskell, 

1953),  follows  Harkrider  (1964)  for  the  solution  of  the  wave 
equation  in  an  n-layered  solid  half  space.  Layer  matrices  of 
the  type  given  by  Dorman  (1962)  for  computing  dispersion  in 
an  n-layered  liquid  -  solid  half-space  are  used  for  the 
liquid  layers.  An  application  of  the  theorem  that  the  inverse  of  the 
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product  matrix  for  the  layered  system  above  the  source  has  the  same 
form  as  the  inverse  of  a  layer  matrix  reduces  the  integrand 
of  the  integral  solution  to  a  simple  form  in  terms  of  product 
matrices  derived  in  the  source-free,  plane-wave  case.  The 
integral  over  wave  number  from  zero  to  infinity  has  singular¬ 
ities  in  the  integrand  and  is,  therefore,  conveniently  trans¬ 
formed  into  the  complex  c  =  k  +  ix  plane.  Performing  the 
contour  integration  in  the  c  plane  and  noting  that  branch 
line  integrals  corresponding  to  branch  points  generated  by 
each  layer  matrix  except  the  last  are  zero,  we  are  left  with 
an  expression  analogous  to  that  of  Sorensen  (1959)  and  Leslie 
and  Sorensen  (1961)  for  the  two-layer  liquid  half-space  with 
a  high-speed  bottom.  By  the  proper  choice  of  contours, 
complex  poles  are  displaced  to  an  unused  sheet  of  the  two-leaved 
Riemann  surface,  and  the  integral  solution  for  the  multilayered 
system  reduces  to  a  summation  of  normal  modes  plus  the  sum  of 
two  integrals  over  wave  number,  one  along  the  real  5  axis  and 
the  other  along  the  imaginary  c  axis.  Both  of  these  integrals 
are  conveniently  computed  by  a  Gaussian  quadrature  formula. 

The  physical  interpretation  of  the  final  solution  is 
straightforward.  The  normal-mode  terms  correspond  to  waves 
trapped  in  the  Arctic  sound  channel;  that  is,  refracted  and 
surface  reflected  (RSR)  sounds  and  reflected  sounds  incident 
on  the  bottom  beyond  the  critical  angle.  The  integral  over 
wave  number  along  the  real  axis  corresponds  to  waves  incident 
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on  the  bottom  at  angles  greater  than  the  critical  angle  of  total 
reflection.  Since  energy  leaves  the  guide  continuously  as  waves 
travel  down  the  guide,  this  term  is  of  importance  only  to  ranges 
of  about  ten  times  the  water  depth.  The  contribution  to  the  sound 
field  of  the  integral  along  the  imaginary  axis  is  of  importance  only- 
very  near  the  source  since  the  integrand  decays  exponentially  with 
range  and  wave  number. 

The  usefulness  of  the  normal  mode  terms  for  describing  long- 
range  explosive  sound  transmission  in  the  central  Arctic  Ocean 
is  shown  by  Kutschale  (1969) .  In  that  work  the  present  formulas 
for  the  normal  modes  were  extended  to  explosive  sources  and  the 
Fourier  integral  for  each  mode  was  evaluated  by  the  principle  of 
stationary  phase  (Pekeris,  1948).  Attenuation  by  the  rough  ice 
boundaries  was  also  included,  although  omitted  here. 

FORMAL  SOLUTION 

Source  Free  Case 

Consider  the  n-layered  interbedded  liquid-solid  half-space 
shown  in  Figure  1  in  which  the  last  layer  of  infinite  thickness  is 
liquid  and  has  a  higher  sound  velocity  than  liquid  layers 
immediately  above  it;  that  is,  a  high-speed  bottom.  A  point 
source  of  harmonic  waves  is  located  in  one  of  the  liquid  layers. 

The  velocity  potential  in  the  m-th  layer,  which  is 
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assumed  to  be  liquid,  satisfies  the  wave  equation 
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Substituting  expressions  for  [Am  -  A^]  and  [Am  +  A^]  from 
equations  (1)  in  equations  (2)  yields 


C  p^vv\ 
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In  general  for  an  (\-layered  liquid  half-space 

%  • 


or 


o 


«  •  » 


ia>\,  uo,; 


since  the  vertical  particle  velocity  and  pressure  are  continuous 
across  each  interface.  For  the  n-th  layer 
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Applying  boundary  conditions  that  the  pressure  vanishes 
at  the  surface  and  that  no  upgoing  waves  travel  from 
infinity,  (PZZ)Q  =  0,  Ar  =  0,  the  relation  (3)  reduces  to 
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Consider  now  a  solid  layer  between  two  liquid  layers  or  at 


the  surface  of  the  laminated  half space.  For  this  layer  we  may 
write  (Thomson,  1950;  Dorman,  1962  )  « 
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where  the  fourth  element  of  the  column  vector  is  the  tangential 
stress  which  is  zero  at  a  solid-liquid  boundary.  The  matrix 
elements  are  given  in  the  appendix  and  may  be  derived  following 
Haskell  (1953).  Equation  (4)  yields  the  three  equations. 
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from  equation  (5)  in  equations  (6)  and  (7) 


yields 


(A<vv7)a 
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which  is  the  same  as  the  matrix  relation  for  a  liquid  layer. 
In  general,  then,  for  interbedded  solids  and  liquids 
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Point  Source  of  Harmonic  Waves  in  a  Liquid  Layer . 

Divide  the  source  layer  into  two  layers  as  shown  in 
Figure  1.  At  z,=D  the  pressure  is  continuous.  The  vertical 
particle  velocity  is  continuous  everywhere  in  the  plane  defined 
by  z  =D  except  at  the  point  source  where  the  liquid  above  and 
below  the  source  moves  in  opposite  directions.  This  may  be 
expressed  by  writing 

-  3  K . 

In  matrix  notation  for  the  source  layer 


For  the  layers  below  z  =D  we  have 


and  for  the  layers  above  z =D 
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Returning  now  to  the  relation  for  the  n-th  liquid  layer 


(11)  may  be  written 
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Harkrider  (1964)  has  shown  that  for  layered  solids  the  inverse 
of  the  product  matrix  has  the  same  form  as  the  inverse  of  the 
layer  matrices.  Employing  an  extension  of  this  theorem  to  inter- 
bedded  liquids  and  solids  we  get 


Hence  (12)  may  be  w rltten  as 
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Adding  equations  (1*0  we  get 


0  =  (J]_]_  +  ^21^  ^  +  ^Ji2  +  ^22^  ^  and  ^  We 


let  T  =  J^]_  +  2 2\  an<^  ^ 


+  and  solve  for  X  we  have 


(15) 


21. 


Carrying  out  the  matrix  multiplication  J  =  E-"1  A,  expressions  for 


T  and  V  are  given  by 
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Hence,  from  the  second  equation  of  (13)  and  equations  (16) 
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Therefore  the  integral  solution  for  the  surface  verticle  particle 
velocity  is  ^  ^  i  ^  k  — 
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Evaluation  of  the  Integral  Solutions. 


The  integral  solutions  have  singularities  corresponding 


to  zeros  of 
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These  are  simple  poles  for  complex  wave  number  and ,  therefore, 
the  integral  solutions  are  conveniently  evaluated  by  contour 
integration  in  the  complex  ^  =  k  +  i  ^  plane.  The  contours 

are  shown  in  Figure  2.  First  the  transformation 
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for  each  solid  layer.  But  since  the  integrands  are  even 
functions  of 
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fkp1  -  o 


considered . 


I  •  I  i  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  1  I  '  I  n  I  '  I  I  I  '  I  '  I  »!■ 


25. 


The  Reimann  surface  h?s  two  sheets.  To  satisfy  the  vanishing 
of  UlT  as  C  approaches  infinity,  we  must  remain  on  the  sheet  of  the 
Riemann  surface  where  the  real  part  of  C.  is  greater  than  zero. 

For  the  contours  shown  in  Figure  2  complex  poles  of 

—^==7  +  -a,v-  ■=  O 

have  been  displaced  to  an  unused  Reimann  sheet  and  all  poles 

are  on  the  real  k— axis.  See  Ewincrf  Jardetzkv  and  Press  (1957  , 
Sages  135  to  137)  for  the  proof  for  a  two-layer 

liquid  half  space.  Along  the  branch  cut 
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From  the  formulas  Hp  c t r)  =  2~  K0C  tP)  *"a 
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where  the  layer  matrices  CX /yv><  are  defined  in  the  Appendix. 
There  are,  as  before,  simple  poles  of  I2  along  the  real  axis  at 


for 


This  is  the  period  equation  for  the  laminated  half  space. 
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In  summary,  the  expressions  for  vertical  particle  velocity 
at  the  surface,  vertical  particle  displacement  at  the  surface, 
and  pressure  in  a  liquid  layer  are: 
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has  been  used. 
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For  programming  it  is  convenient  to  write: 


ur0  -  -Y/^V*  -  c V/3+  c V/^  -  l 


and  for  the  absolute 


value  of  vertical  particle  velocity,  a  quantity  conveniently 

measured  in  transmission  experiments, 
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Likewise 


for  vertical  particle  displacement 
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and  for  pressure 
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The  layer  matrices  for  the  integral  along  the  real  axis 
are  programmed  as 


for  solid  layers  and 

(a.  ^ 

-  w 

for  liguid  layers.  For  the  integral  along  the  imaginary  axis 
the  matrices  are 


for  liquid  layers.  Likewise,  for  the  ndrmal  mode  terms  the 
layer  matrices  are  programmed  as 


39. 


for  solid  layers  and 


for  liquid  layers. 
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NUMERICAL  COMPUTATIONS 


Computer  programs  were  written  in  double  precision  Fortran 

2 

IV  to  compute  the  pressure  in  dynes/cm  and  the  vertical  particle 
displacement  in  millimicrons.  The  programs  are  run  on  the  IBM 
360/91.  For  numerical  results  presented  here,  computing  time 
was  under  six  minutes.  It  appears  that  practical  limits  of 
numerical  precision  make  the  present  development  most  useful  for 
the  Arctic  at  frequencies  below  50  Hz  in  water  depths  to  1  km. 

At  higher  frequencies  or  in  greater  water  depths  the  integrands 
may  be  so  oscillatory  that  it  is  difficult  to  achieve  the  desired 
accuracy  in  the  numerical  integrations  without  excessive  computing 
time . 

Computations  are  made  in  two  stages.  The  first  program,  an 
extension  of  Dorman's  (1962)  PV  7  dispersion  program,  computes 
phase-  and  group-velocity  dispersion,  the  excitation  function  depen¬ 
dent  only  on  layering  of  the  medium,  and  the  excitation  function  for 
the  particular  source  and  detector  depths.  In  the  first  case  for 
the  m-th  normal  mode  the  excitation  function  is  defined  by 
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for  vertical  particle  velocity.  In  the  second  case  the  excitation 
function  is 
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for  pressure  and 
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for  vertical  particle  velocity.  These  definitions  were  chosen  to  be 

ful  also  at  long  ranges  where  VC  Ik  c)  = 

(% 


use 
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va-kr) 


The  second  program  computes  the 
three  integrals,  the  normal  mode  contribution  to  the  sound  field, 
and  writes  and  plots  the  absolute  value  of 
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\ 


,  or  wQ  as  a 
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TABLE  1. 


Model  Parameters  for  Sorensen's  (1959)  Computations:  Source 

2 

Pressure  Amplitude  re  lm,.3048  dyne/cm  ;  Source  Depth,  15.2  m; 
Hydrophone  Depth,  15.2m;  Range,  3.048  m 


Layer  Thickness,  m  Compressional 

Velocity,  m/sec 


Shear 

Velocity,  m/sec 


Density 
gm/ cm0 


18.3 


1463.0 


0 


1.03 


Infinite 


1609.3 


0 


1.24 
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TABLE  2. 


Integral  Real  Axis,  Real  Part 


F  requency, 
Hz 


Sorensen’s  (1959)  Value, 
dynes/cm2 


Our  Value, 
dynes/cm2 


10 

.004713 

.004685 

20 

.030170 

.030153 

40 

.040278 

.040485 

80 

.023193 

.023121 

160 

-  .021839 

-  .022201 

320 

-  .074395 

-  .074011 
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TABLE  3. 

Integral  Real  Axis,  Imaginary  Part 


Frequency, 

Hz 


Sorensen1 s  Value, 

dynes/cm^ 


Our  Value, 

dynes/cm^ 


10 

-  .002567 

-  .002551 

20 

-  .022260 

-  .022240 

40 

-  .061960 

-  .062274 

80 

-  .050986 

-  .050641 

160 

-  .047979 

-  .047750 

320 

.006299 

.008316 
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TABLE  4. 


Integral  Imaginary  Axis 


Frequency, 

Hz 

Sorensen’s  Value, 
dynes/cm2 

Our  Value 

dynes/cm2 

10 

.095080 

.095146 

20 

.081938 

.082076 

40 

.041887 

.041759 

80 

.030795 

.030937 

160 

.030214 

.030356 

320 

.010593 

.011043 

DEPTH  (M) 


H6. 


SOUND  VELOCITY  (M/SEC) 


Fig.  3 


TABLE  5. 
Model  A 


Layer 

Thickness,  m 

Compressional 
Velocity,  m/sec 

Shear  Velocity, 
m/sec 

Density, 

gm/cm3 

50 

1431).  0 

0 

1.03 

50 

1437.0 

0 

1.03 

50 

1440.0 

0 

1.03 

50 

1443.0 

0 

1.03 

50 

1446.0 

0 

1.03 

50 

1450.0 

0 

1.03 

50 

1453.2 

0 

1.03 

Infinite 

1600.0 

0 

1.20 
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function  of  range.  In  this  program  provision  is  also  made  to 
compute  the  three  integrals  as  a  function  of  frequency  or  detector 
depth  at  a  specified  range.  The  integrals  are  evaluated  by 
Gaussian  quadrature  formulas  (see,  for  example,  Hildebrand,  1956). 
Care  must  be  taken  to  employ  a  sufficient  number  of  Gaussian  points 
to  obtain  the  desired  accuracy  of  the  integrals.  A  thirty-two 
point  Gaussian  quadrature  formula  was  found  to  give  sufficient 
accuracy  for  the  computations  presented  here  out  to  a  range  of 
1500  m,  but  at  longer  ranges,  higher  frequencies,  or  greater  water 
depths  a  ninety-six  point  formula  is  used.  Abscissas  and  weight 
factors  for  the  Gaussian  integration  are  given  by  Davis  and 
Rabinowitz  (1956,  1958).  The  integration  along  the  imaginary 
axis  extends  to  infinity,  but  in  practice  a  finite  upper  limit  is 
chosen  which  gives  sufficient  accuracy.  This  may  be  done  because 
the  integrand  is  an  exponential  function  of  -kr  for  large  k  or  r. 

The  programs  were  checked  by  computations  for  the  same  two-layer 
model  of  Table  1  used  by  Sorensen  (1959).  Tables  2,  3  and  4  show 
that  the  two  sets  of  numerical  results  are  in  very  close  agreement. 

As  an  illustration  of  the  progrms  for  the  Arctic,  computations 
were  made  for  the  layered  model  A  of  Figure  3 ,  which  closely  follows 
the  observed  variation  of  sound  velocity  with  depth.  Additional 
parameters  for  the  model  are  given  in  Table  5.  This  is  a  prelim¬ 
inary  model  and  the  ice  sheet,  which  averages  about  3-m  in  thick¬ 
ness  in  the  central  Arctic  Ocean,  was  not  included.  The  effects  of 
the  ice  sheet  on  propagation  are  currently  under  investigation. 
Figure  A  shows  phase  -  and  group-velocity  dispersion  and  the 
excitation  function  for  pressure.  At  a  frequency  of  10  Hz  only 
two  normal  modes  are  excited.  The  range  dependence  of  the  absolute 
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value  of  pressure  for  the  normal  mode  contribution,  the  integral 
contribution,  and  the  total  pressure  are  shown  in  Figures  5  through 
10. 
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APPENDIX 

Matrix  elements  for  normal  modes  (see  Haskell,  1953,  and 
Dorman,  1962). 

Solid  layers 

=(5^,,= cos,^ 
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Matrix  elements  for  integrals  along  real  axis: 
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Liquid  layers 
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Matrix  elements  for  integral  along  imaginary 
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Liquid  layers 
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Corrections 


Page  7,  line  3,  should  read  as  follows: 
in  terms  of  elements  of 

Page  13,  line  4,  should  read  as  follows: 


Page  16,  line  8,  should  read  as  follows: 


■&n-,  ^ 

! 

&<T\ 

■  b^a^-'-b, 

o 

Page  19,  line  1,  should  read  as  follows: 

A  ~  Ac,^  ^n-\^n-a  '  ’  '  '°n-\o 

Page  22,  line  7,  should  read  as  follows: 

at  depth  #  at  the  bottom  of  the  layer  is 


Page  22,  line  10,  should  read  as  follows: 

<Xt>  -  \  •  •  •  Q  t 

Page  28,  line  2,  should  read  as  follows:  but  from  the  relation 

A/\  '-I  )  ^V'  ~  ^'3.  Aq,  “  l 
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Page  41,  second  line,  should  re<  as  f ollows : 
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